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We numerically study the temporal evolution of a black-hole laser configuration displaying a pair 
of black and white hole horizons in a flowing atomic condensate. This configuration is initially 
prepared starting from a homogeneous flow via a suitable space-dependent change of the interaction 
constant and the evolution is then followed up to long times. Depending on the values of the system 
parameters, the system typically either converges to the lowest energy solution by evaporating away 
the horizons or displays a continuous and periodic coherent emission of solitons. By making a 
physical comparison with optical laser devices, we identify the latter regime of continuous emission 
of solitons as the proper black-hole laser effect. We include some movies of the temporal evolution 
of the spatial density and velocity profiles in the most significant cases. 

PACS numbers: 03.75.Kk 04.62.-|-v 04.70.Dy 


I. INTRODUCTION 


Since the pioneering work by Corley and Jacobson [ij, 
gravitational configurations showing a pair of neighbor¬ 
ing horizons have attracted a great deal of attention: for 
a quantum Bose field with a superluminal dispersion, the 
negative energy partner of the Hawking emission by the 
outer horizon can bounce on the inner horizon and travel 
back to the outer one, so to stimulate further Hawk¬ 
ing emission. In suitable conditions, this process may 
give rise to a dynamical instability, the so-called black- 
hole laser instability, and then to the exponential growth 
of a self-amplifying coherent Hawking emission. Most 
remarkably, this coherent emission is expected to have 
very different statistical properties from the spontaneous 
Hawking radiation from a single horizon, which originates 
from the parametric amplification of zero-point quantum 
fluctuations and whose two-mode squeezed state nature 
reduces to a simple thermal state upon tracing out the 
Hawking partner emission [H- 

While this intriguing emission seems to have little 
chances to be observed in astrophysics, it is central to 
present-day studies of optical and condensed matter ana¬ 
logues of gravitational systems. Such analog models, first 
proposed by Unruh Q as table-top systems where to 
study the analog Hawking emission process, are presently 
attracting a great interest from both the theoretical and 
experimental points of view. Among the specific sys¬ 
tems under most active study, we can mention nonlinear 
optical systems [3], surface waves in water tanks ii, 
quantum fluids of light and Bose-Einstein con¬ 
densates of ultracold atoms A clear evidence of 

the presence of a black-hole horizon has been obtained in 
both fluids of light m and atomic gases [l^ and strong 


experimental efforts are presently being devoted to the 
detection of the spontaneous Hawking emission. A re¬ 
markable first experimental evidence of a black-hole laser 
instability has been recently reported in [l^ by looking at 
the fast growth in time of a complex density modulation 
pattern in between the black- and white-hole horizons. 

Whereas the first works on black-hole laser configu- 
rations in astrophysics [I| and in analog models [l^ 
0 have focussed on the linear dynamics at early times 
after the onset of the instability, the most interesting 
physics occurs at late times when the exponentially grow¬ 
ing black-hole laser emission has become strong enough 
to exert a measurable back-reaction effect onto the hori¬ 
zons. Focusing on flowing atomic condensates for which 
theoretically accessible multi-horizon configurations were 
proposed along the lines of [H, pioneering recent 
works [2l|, [ 2 ^ have identified some most remarkable be¬ 
haviors that the system can take in response to the back- 
reaction effect. The black-hole laser configuration can 
quickly relax towards a horizonless sub-sonic flow by 
“evaporating” away the horizons, and, in strongly unsta¬ 
ble configurations, spontaneous oscillations can appear 
and eventually result in the periodic emission of solitons. 

In this Article we extend those pioneering works with 
an extensive campaign of numerical simulations of the 
long-time dynamics of atomic condensates starting from 
a black-hole laser configuration with a pair of neigh¬ 
boring black- and white-hole horizons. Our simula¬ 
tions are based on the numerical integration of the one¬ 
dimensional time-dependent Gross-Pitaevskii equation 
describing the condensate dynamics at mean-field level. 
Besides confirming the qualitative behaviors anticipated 
in Refs, [l^, a thorough classification of the different 
regimes that the system can display depending on the 
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system parameters is presented and some physical expla¬ 
nation is given for the more complex and intriguing ones. 

A special attention is devoted to the spontaneously os¬ 
cillating regimes where the classical dynamics of the sys¬ 
tem quickly forgets its initial condition and approaches 
a limit cycle attractor at late times; as a consequence 
of these oscillations, solitonic sound waves keep being 
emitted into the surrounding condensate for indefinite 
times in a continuous and periodic way. Even though the 
“black-hole laser” expression has been used in the past to 
refer to different aspects of this physics, a detailed com¬ 
parison to continuous-wave optical laser devices suggests 
that this self-oscillating regime is the closest counterpart 
of laser operation in quantum optics, where the periodic 
oscillation of the electromagnetic field in a laser cavi^ 
leads to the emission of monochromatic radiation [2l-Q. 
The clarification of the relation between the laser oper¬ 
ation in quantum optics and the black-hole laser phe¬ 
nomenon in gravitation and analog models is another 
main contribution of this work. 

The scheme of the paper is the following. In Sec. [II] 
we revisit the basic theory of black-hole lasers in Bose- 
Einstein condensates. In Sec. uni we present our numeri¬ 
cal results and we classify the different regimes according 
to the number of linearly unstable modes in the initial 
condition (Secs. IIII A H and Sec. IIIIA2I) . The continuous 
soliton emission is then discussed in Sec. IIII Bl Compar¬ 
ison between this mechanism of emission of solitons with 
the operation of optical laser devices is given in Sec. IIVI 
Conclusions and outlook are finally given in SeclV) Tech¬ 
nical details about the method used for the numerical 
calculations are presented in Appendix!^ A summary of 
the movies and of the parameters used in each of them 
is given in Appendix [B] 


II. BLACK-HOLE LASER CONFIGURATIONS 

In this section, we briefly revise the basic concepts of 
the theory of black-hole lasers in flowing Bose-Einstein 
condensates showing two neighboring horizons of oppo¬ 
site black- and white-hole nature. These concepts are 
instrumental to physically understand the results of the 
numerical simulations that will be discussed in the fol¬ 
lowing sections. For further details, we refer the reader 
to the recent literature on the topic, in particular Refs. 

[TMl,!!!!!!. 


A. The system and the theoretical model 


We consider a flowing atomic one-dimensional Bose- 
Einstein condensat e 1^ n ear T = 0 in the so-called ID 
mean-field regime [2J, l25l | where the transverse confine¬ 
ment length is much larger than the s-wave atom-atom 
scattering length and the transverse wavefunction is not 
distorted from its Gaussian shape in the absence of in¬ 
teractions. Under these assumptions, the condensate 


dynamics is accurately described in terms of a macro¬ 
scopic wave function which evolves according to a time- 
dependent one-dimensional Gross-Pitaevskii (GP) equa¬ 
tion: 


ih 


~dt 




+ , ( 1 ) 


where we have allowed for a spatial dependence of the ex¬ 
ternal potential V{x) and of the effective one-dimensional 
atom-atom interaction constant g{x). In experiments, 
both quantities can be controlled and manipulated using 
standard tools of atomic physics [2^. 



FIG. 1; Scheme of the theoretical model here considered for 
studying analog black-hole laser configurations, where the 
stationary GP wave function is a plane wave of the form 
5'(a:) = ydioe*'^^. Upper panel: Scheme of the spatial profile 
of the constant coupling g{x) (solid blue) and the external po¬ 
tential V{x) (dashed red), which are piecewise functions sat¬ 
isfying Eq. Lower panel: Scheme of the velocity profile. 
The local speed of sound, c(x) = ^g{x)no/m, (solid blue line) 
crosses twice the uniform flow velocity, v{x) = v = hq/m > 0, 
directed in the rightwards direction (red dashed line). The 
size of the internal supersonic region between the two hori¬ 
zons is X and the supersonic sound speed is C 2 < v. 

Among the different configurations proposed and/or 
experimentally realized to implement analog models of 
gravity in flowing atomic condensates, a theoretically 
very convenient one (although very idealized from the 
experimental point of view) was introduced in [3113 
and then widely used for analytical and numerical stud¬ 
ies of black holes and black-hole lasers [3.13. 13. [3| . 
The idea is to have piecewise constant V and g satisfying 
at all points the condition 

g(x)no+ V(x) = Eb , (2) 

Eh being a constant. Under this condition, a solution of 
the GP equation exists with the condensate keeping at 
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all times a plane wave shape 

(3) 

with a homogeneous density no, a speed v = hq/m, and 
a chemical potential ^ = Et + li^q^/2m. 

As illustrated in Fig. [1] we adopt a piecewise homo¬ 
geneous configuration with g{x) = gi and V(x) = Vi for 
|a:| > X/2 and g{x) = 52 and V{x) = V 2 for |a:| < X/2. 
The spatial variation of g{x) results in a corresponding 
variation of the sound speed, 


lasing. In the following subsection we shall present the 
structure of the different nonlinear stationary states of 
the GP equation for the given spatial profiles of g{x) and 
V(x). These states play a central role in determining the 
long time behavior of our system. 

To study the dynamical stability of the initial plane 
wave stationary configuration, one can consider linear 
perturbations of the wave-function around the station¬ 
ary plane wave solution, 

t) = [1 + t)] (5) 


c(a;) = Cl,2 = 


5i,2«o 


( 4 ) 


for respectively |x| ^ X/2. A black-hole laser configu¬ 
ration appears when the external regions are subsonic 
and the internal region of size X is supersonic, i.e., 
C2 < V < Cl. In this case, two sonic horizons are in fact 
formed, the upstream one being located at a: = —X/2 
with a black-hole nature, the downstream one being lo¬ 
cated dX X = X/2 with a white-hole nature. The match¬ 
ing condition ([2]) imposes that Vi -f giu-o = V 2 + 

All our numerical simulations of the GP equation start 
from the spatially homogeneous plane wave state ©: 
in an actual experiment, such an initial condition could 
be implemented by applying a spatially-selective sudden 
jump of the interaction constant and of the external po¬ 
tential in the inner region |a;| < X/2 of an initially ho¬ 
mogeneously flowing condensate. A similar configuration 
has been adopted in recent numerical studies of Hawk¬ 
ing radiation from black hole configurations [l^. Even 
though this configuration might be extremely challeng¬ 
ing to realize in the lab, it allows a simpler theoretical 
understanding of the underlying physics easier since it 
eliminates unwanted emission channels as condition @ 
guarantees that the horizons do not generate any deter¬ 
ministic excitations of the condensate, e.g. Bogoliubov- 
Cherenkov waves (ttI . . For the same reason, we will 
also assume that the condensate is infinitely extended on 
either side of the supersonic region. 

For computational convenience, hereafter we rescale 
the wave-function by the homogeneous density no, 
'I'(a;,t) —>■ and adopt units such that h = 

TO = Cl. Within this convention, the initial homogeneous 
condensate density |'I'(a;,t = 0)p and the sound speed ci 
are 1, the flow velocity is v = q and hence the system is 
completely characterized by the three values of (c 2 , c, X). 


B. Unstable modes 


To introduce the most relevant features of the system 
investigated in this Article and to facilitate the interpre¬ 
tation of the outcome of our numerical simulations, we 
will give a quick summary of the results of Refs. (^.1^. 
In this subsection, we briefly describe the structure of 
the (linear) dynamical instabilities governing the tran¬ 
sient dynamics during the initial stages of the black-hole 


whose dynamics is described by the Bogoliubov-de 
Gennes (BdG) equations [ 2 ^. In our piecewise homo¬ 
geneous geometry, the solutions of the BdG equations 
for the perturbation (5'I'(a;,t) can be written within each 
subsonic or supersonic region as linear combinations of 
plane waves satisfying in the laboratory frame the Bo- 
goliubov dispersion relation 

^4 

(w — vk/^ = Cj -I- —, z = 1, 2 (6) 

For a given frequency w, the global solutions to the BdG 
equations are obtained by matching the plane-wave solu¬ 
tions at the boundaries between the external and internal 
regions. 

While single-horizon configurations of both black- and 
white-hole nature reduce to a scattering problem on a sta¬ 
ble horizon [l^, [ 2 ^, [ 2 ^ , the physics of two-horizon config¬ 
urations is much more intriguing [Hill,[mill,[13: if the 
size X of the supersonic region is sufficiently large, dy¬ 
namical instabilities occur, signalled by the appearance 
of solutions with a complex frequency 7 „ = -I- zF^ 

with a growth rate r„ > 0. The appearance of such dy¬ 
namically unstable modes provides the black-hole laser 
effect. 

As X increases further, more and more unstable modes 
appear, labeled by the quantum number n = 0,1, 2.... In 
particular, the critical length at which a new dynamical 
instability appears for given (w,C 2 ) is [2l|: 

Xn = Xq + uXq (7) 


where 


arctan 


Xn = 




and Ao = 




( 8 ) 


The intermediate values X = Ar„_|_i/2 separate qualita¬ 
tively different unstable behaviors: while for Xn < X < 
Xn+ 1/2 the complex frequency 7 ^ is purely imaginary, for 
-^n-i- 1/2 < X < Xn+i it acquires a finite oscillation fre¬ 
quency ujn yf 0. Remarkably, the instability rate drop s to 
r„ = 0 right at the transition point X = Xn^ii 2 


C. Non-linear stationary solutions 

While the short-time dynamics of the system is well 
captured by the BdG equations, the linear approximation 
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underlying them eventually breaks down at some point 
after the onset of the instability and one needs to go 
back to the full time-dependent GP equation ([T]). In 
particular, a step of crucial importance to understand 
the long-time behavior of the system is the identification 
of the different stationary states of the system. 

We review now the inhomogeneous non-linear GP sta¬ 
tionary solutions that exist for the same piecewise profile 
of g{x),V{x) considered above. Following Ref. [^, we 
specifically focus on those solutions that asymptotically 
match the homogeneous plane wave solution, 'I'(a;,t) = 
The quest of such particular set of solutions 
is motivated by the fact that the total energy and num¬ 
ber of particles are conserved during the growth of the 
instabilities and, thus, only these solutions present a fi¬ 
nite energy and particle number difference with respect 
to the initial plane wave. Hence, these solutions are good 
candidates as ending points of the time evolution of the 
black-hole laser. 

Writing the stationary GP solution in the form 

■^{x,t) = with 4'(a;) = (9) 


given wave function 4' is: 


n[^] = noJdx + (Vix) - 

(13) 



X 


gives us two equations for respectively the phase and the 
amplitude: 


J 

g.A{x) 


A‘^{x)(j)'{x) 

A''{x) J2 
2 2A3{x) ^ 

V{x)A{x) + g(x)A^{x) 


( 10 ) 


( 11 ) 


where the conserved current J is constant and the prime 
symbol ' denotes spatial derivative. For inhomogeneous 
solutions, we can define local flow and sound velocities, 
given by v{x) = (j)'{x) and c(x) = ^^/g{x)p[x), with 
p{x) = A^{x) the density. As the solution is asymp¬ 
totically a plane wave, one has J = v. 

Within each (subsonic or supersonic) region, the cou¬ 
pling constant and the potential are homogeneous and 
we can integrate (I10mi|) to obtain: 


„ C-* 

~ + “Y 


( 12 ) 


where Ci is a conserved quantity for each z = 1,2 region 
and Pi = p—Vi, with Vi the constant value of the external 
potential in the region i. In our units, gi = cf and pi = 
Qi + see Eq. H. In particular, we have gi = 1, 

Pi = l + z;V2 and Ci = l + 2v^. 

As a criterion to choose among the many possible 
solutions of these equations, it is reasonable to expect 
that the initially dynamically unstable system will evolve 
towards a new stationary solution which minimizes its 
grand-canonical energy defined as = H — pN in terms 
of the the Hamiltonian H and the total number of par¬ 
ticles N. The expression for the functional H[4'] for a 


FIG. 2: Density profile of the n = 0, 1,2 nonlinear station¬ 
ary solutions of the GP equation, Eq. (O, for parameters 
V = 0.9, C 2 = 0.5 and X = 10. Dashed vertical black lines 
represent the limits of the internal supersonic region. 


When 41 is a stationary solution of the form ([S]), it is 
an extreme of the previous functional. In particular, the 
energy difference of a stationary solution with respect to 
the initial homogeneous plane wave solution (with energy 
ilh) is 


AH = H[^] -nh = -no y dx ^ [A^ix) - l] (14) 

As the density p{x) is asymptotically fixed in the external 
subsonic regions, an increase of the density within the 
internal |x| < X/2 region leads to an overall reduction of 
the grand-canonical energy. This increase of the density 
carries a finite increment of the number of particles: 


AN = 


dx [p[x) — 1] 


(15) 


In particular, the family of stationary solutions with low¬ 
est grand-canonical energy AH < 0 corresponds to the so- 
called sh-sh solutions . They consist of a shadow- 

soliton solution in the external regions and of elliptic 
functions in the internal region. Specifically, the station¬ 
ary GP wave function 4'„(a:) of these solutions is given 
by: 
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FIG. 3: Sound (solid blue) and flow velocity (dashed red) profiles corresponding to the n = 0,1, 2 stationary GP solutions of 

Fig.H 


^1 - coth 

vl/„(x) = vl/j^"‘(x), |x| < X/2 




y/l — v'^{Ax —-— x) 


X 


(16) 


«'„(a;) = ( y/l - coth 

coth“^ 




\/l — u^(Ax —— + x) 


X 

-iv\ , x> Y 


Ax = 


1 — 1?2 


\/r^ 


= p{-X/2) = p{X/2) 


with the phase Ac/) chosen such the wave function is continuous at a; = ±A/2. The wave function in the internal 
region, 'F“‘(a;), reads 


Pn\x) = ri + (r 2 - ri)sn^(c 2 \/r 3 - nx + (n + l)K{v),v), v = 


(17) 


^T{x) = 


riC2y/r3 - ri 


r2 - ri 
rs - ri ' 

r2 


ri 


n ( am(c 2 \/r 3 — rix + (n + l)K{i'), v), 1-z/ ) — 11 f am ((n + l)K{v),v) , 1- -,v 


ri 


In the above equations, K{v) is the complete elliptic inte¬ 
gral of the first kind, n((/), v) is the incomplete elliptic 
integral of the third kind and sn is the corresponding 
Jacobi elliptic function [s^, [Slj. The quantum number 
n = 0,1,2... characterizes the solution and gives the 
number of density minima inside the supersonic region. 
The parameters 0 < ri < r 2 < rg are the roots of the 
following polynomial equation for the density: 

g2P^ - 2/i2p^ -h C 2 P - = 0 (18) 

The specific value of C 2 is obtained by imposing the conti¬ 
nuity of the amplitude of wave function and of its deriva¬ 
tive at a; = ±A/2. On one hand, after using Eq. (|T^ . 
we get: 

Pm = ^ + (19) 


On the other hand, using Eq. (HU), one finds 


A = 


2 

C2\/r’3 - ri 


(n + l)K{v) — sn ^ 



( 20 ) 

This equation gives an implicit equation for C 2 as pm 
and the values of the three roots ci 2,3 depend on C 2 . 
The minimum length at which a new solution appears 
is found in the limit —>■ 0, where the Jacobi elliptic 
functions reduce to the usual trigonometric functions, 
sn(it, z^) i:: sin(rz). Using this relation, it can be easily 
shown that the minimum length at which the n non¬ 
linear solution appears coincides with the onset of the 
dynamical instability at A = A„, as given by Eq. 0: 
there is therefore a remarkable correspondence between 
dynamically unstable modes and non-linear solutions of 
the GP equation, as shown in Ref. [2l|. 


The density profile p{x) = A^{x) of the first few non¬ 
linear solutions of Eq. dSD for a given configuration is 
plotted in Fig. HI The grand-canonical energies of these 
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solutions satisfy Af2n< Afii < Ari 2 ... < 0. As first 
anticipated in Refs. j2l|, by studying linearized fluc¬ 
tuations around the stationary state, only the n = 0 
solution is dynamically stable while all other solutions 
are dynamically unstable: a time-dependent integration 
of the GP equation fully confirms this fact, see the next 
Section. 

A qualitative argument justifying the fact that the n = 
0 solution is the most stable one can be put forward as 
follows: the density profile of this solution (displayed in 
solid black line in the left panel of Fig. la shows an 
accumulation of atoms in the central region at |a;| < X/2, 
so that the local sound velocity increases and, according 
to mass conservation at a fixed current, the flow velocity 
is reduced. As a result, the supersonic unstable character 
of the central region is reduced and, for sufficiently large 
values of A, is suppressed so the flow becomes everywhere 
subsonic, restoring in this way the full dynamical stability 
of the configuration. We can observe this behavior in the 
left plot of Fig. [31 where the sound and flow velocity 
profiles for the n = 0 solution of the left plot are depicted. 

The situation is of course completely different for the 
higher n > 0 nonlinear solutions shown in the central 
and right panels of the same figure: the density min¬ 
ima that are present in these solutions in the central 
(|a;| < X/2) region always correspond to a locally su¬ 
personic flow; in particular, for this specific configura¬ 
tion, the n = 2 solution represented in the rightmost 
panel is completely supersonic in the whole central re¬ 
gion. As usual in hydrodynamics the presence 

of such localized super-sonic flow regions in combination 
with spatial inhomogeneities is typically responsible for 
dynamical instabilities of the flow pl2H3^ . 


III. NUMERICAL RESULTS 


After the brief review of the general concepts of black- 
hole laser physics, we are now able to start presenting our 
numerical results for the time evolution of the black-hole 
laser. As mentioned in the previous section, we consider 
a system with a homogeneous potential and coupling con¬ 
stant and then, at t = 0, we introduce a sudden spatially- 
dependent quench of the coupling constant g(x) and the 
potential V{x) in such a way that they still satisfy Eq. 
©• The corresponding time evolution of the wave func¬ 
tion for t > 0 is governed by the following time-dependent 
GP equation: 


^~dt 


( 21 ) 


where the constant term Eh of Eq. @ has been sub¬ 
tracted from the usual time-dependent GP equation © 
for numerical convenience. We recall that, in our rescaled 
units. Eh = g{x) V(x). In particular, we will concen¬ 
trate on the configuration sketched in Eig. |TJ where g{x) 
and V (x) are piecewise constant functions as described 
around Eq. and the initial condition is taken in the 


form of the stationary plane wave (O, On one hand, 
these choices allow us to make direct use of the general re¬ 
sults reviewed in Secs. Ill Bl and III Cl On the other hand, 
we also expect that most of the main conclusions of our 
work should be transferrable to the experimental config¬ 
uration of [l^, at least for times shorter than the actual 
transit time of the incident condensate on the potential 
barrier. 

To ensure that all the dynamical instabilities are duly 
triggered in the simulations, a weak random noise is 
added on top of the plane wave initial wavefunction. The 
size of the numerical grid is taken sufficiently large to 
avoid finite size effects, Lg ^ 10^ X. Undesired re¬ 
flections from the boundaries of the integration box are 
further suppressed by implementing a diffusive term at 
the edges of the grid. This allows us to extend the nu¬ 
merical simulations up to long times t ~ 10'^ without suf¬ 
fering from numerical artifacts. In the next subsections, 
the main features of the time evolution of the system will 
be characterized as a function of the relevant parameters 
{c 2 ,v,X). A more detailed description of the numeri¬ 
cal integration scheme and of the implementation of the 
diffusive term is given in Appendix [^ 


A. Long-time stationary state 


As a first step, three main regimes can be identified 
depending on whether X < Xq(v,C 2 ), Xq(v,C 2 ) < X < 
Xi{v,C 2 ) or Ai(v, C 2 ) < X with Aq.i defined in ©■ 

The first case X < Xo{v, C 2 ) is trivial as there are no 
instabilities and the time evolution reduces to the dy¬ 
namics of the (very weak) noise on top of the stationary 
plane wave solution. In the following we therefore focus 
our attention on the other two cases that show the most 
interesting physics. 

Before entering into the details, some preliminary re¬ 
marks are however in order. Eirst, we note that the limit 
C 2 —>■ 0 is ill-defined, as in this limit the n = 0 non-linear 
solution presents an infinite accumulation of particles be¬ 
tween the two horizons, which cannot be achieved start¬ 
ing from a finite condensate as that considered in our 
numerical simulations. We then restrict our simulations 
to finite values of C 2 for which the particle accumulation 
AN is much smaller than the total number of particles 
in the system, AN <C N. 

As a further tool to facilitate the reader in following 
our discussion, we supplement the figures presented in 
the main text of the Article with a few Movies that high¬ 
light the main qualitative features of the black-hole laser 
dynamics. Links to the Movies are included in the main 
text of the Article, while the system parameters chosen 
for each of them are summarized in Appendix |B] 
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1. Single unstable mode Xq < X < Xi 

In the window Xo(n, C 2 ) < X < Xi(n, C 2 ), only one 
unstable mode is present, associated with the appearance 
of the lowest energy n = 0 nonlinear solution described 
by (dll), whose energy lies below that of the homogeneous 
solution. After some transient governed by the unstable 
mode, we can therefore expect that the system will expel 
the extra energy in the form of waves propagating to 
X —>■ ±00 and eventually relax towards this non-linear 
stationary solution. 

Whilst this behavior was indeed observed in Ref. [^ . 
we have noticed that it is restricted to a certain region 
in parameter space, corresponding to sufficiently low val¬ 
ues of V and high values of C 2 . Outside this region, the 
instability of the supersonic region becomes too severe 
and the convergence towards the n = 0 stationary solu¬ 
tion is replaced by a time-dependent regime of continuous 
emission of solitons (CES). Importantly, we have numeri¬ 
cally observed that the choice between the two behaviors 
only depends on the system parameters {c 2 ,v, X) and not 
on the specific initial configuration of the noise used to 
trigger the instability. In this subsection, we focus our 
attention on the first behavior while the detailed char¬ 
acterization of the CES regime and its phase diagram is 
postponed to Sec. IIIIBI 

Eor given values of {v,C 2 ) outside the CES window, 
two cases Xo{v,C 2 ) < X < A'i/ 2 (n,C 2 ) and Xi/ 2 {v,C 2 ) < 
X < Xi{v,C 2 ) can be further distinguished as briefly 
mentioned at the end of Sec lII Bl In the first case, 
Xq(v,C 2 ) < X < Xi^ 2 iv,C 2 )j the frequency of the unsta¬ 
ble mode is purely imaginary and the mode does not os¬ 
cillate. As a result, the density shows for times t > I/Eq 
a monotonic exponential evolution of the form: 

p{x, t) Ki 1 + aoSpo{x)e^°* (22) 

where Spq{x) is the linear density perturbation corre¬ 
sponding to the unstable mode and oq its initial ampli¬ 
tude. Depending on the sign of the real number ag, the 
density will either grow or decrease: as the exponential 
evolution is triggered by a random noise on top of the 
stationary solution, the probabilities of the two instances 
are both equal to 1/2. This fact has been numerically 
checked. 

When the density initially increases, the system 
smoothly reaches the n = 0 ground state solution. The 
increase of the density in the central, initially supersonic 
region is associated with a small emission of waves and 
a small soliton to the upstream region (a: —> — 00 ) in or¬ 
der to conserve the total number of particles N. We can 
observe this behavior in the upper row of Eig. |4l which 
shows the result of a simulation for a choice of param¬ 
eters V = 0.75, C 2 = 0.3 and X = 2 that falls in the 
Xq(v,C 2 ) < X < Xi/ 2 {v,C 2 ) window. The complete time 
evolution can be observed in Movie 1, 

On the other hand, when the density initially de¬ 
creases, the system has to emit a larger soliton to the 
downstream region [x —>■ 00 ) in order to compensate the 


initial decrease in the particle density. Once the soliton 
has been emitted, the system is again free to evolve to 
the n = 0 solution by locally increasing the density in 
the central region. This scenario is shown in the central 
and lower row of Fig. |4l which corresponds to a simula¬ 
tion generated with the same parameters {c2,v, X) but a 
different configuration of initial noise; the complete evo¬ 
lution can be observed in Movie 2, 

For Xi/2{v,C2) < X < Xi{v,C2), the system also 
evolves towards the n = 0 stationary solution. The only 
difference is the transient: as the instability keeps oscil¬ 
lating while growing, the dependence on the initial noise 
condition is no longer relevant and no distinction between 
initially increasing and decreasing cases can be made any 
longer. This behavior is shown in Fig. [SJ see also Movie 3 
for the complete picture. 

Our numerical results then confirm that in this stable 
regime the (very different) transient dynamics does not 
play any significant role on the long time behavior of 
the system, which is only determined by the values of 
{c 2 ,v,X) and in all cases eventually tends to the ground- 
state n = 0 solution. 


2. Several unstable modes Xi < X 

For larger sizes X of the central region, several unstable 
modes are present and, because of that, the dynamics is 
potentially much more complicated. From Eq. © , we 
see that the minimum length at which a new unstable 
mode appears is given by: 

= Xn(v = 1, C 2 = 0) = rnr (23) 

Also in this case, two different behaviors are observed de¬ 
pending on the degree of instability: for sufficiently slow 
speeds v and high values of C 2 , the long-time evolution af¬ 
ter a (sometimes complex) transient tends to a stationary 
state. On the other hand, for sufficiently high speeds v 
and low values of C 2 , a CES regime appears. In this sub¬ 
section we restrict ourselves to the former case, leaving 
the discussion of the CES regime for the next subsection. 

It was shown in Ref. [2l| that the unstable mode with 
largest value of n is typically the dominant one in the 
early evolution given its fastest growth rate r„. This fea¬ 
ture has been qualitatively confirmed in our simulations 
by looking at the spatial profile of the density modulation 
at intermediate times, at its oscillating/non-oscillating 
temporal dependence, and at the instability growth rate. 

For what concerns the later dynamics, the situation 
is quite similar to the one discussed in the previous sec¬ 
tion, the main difference being that several stationary 
solutions with Afl < 0 are now available. While on the 
long run the system typically tries to relax to the lowest 
energy n = 0 solution (which is the only dynamically sta¬ 
ble one) by emitting the extra energy and particles in the 
form of solitons and small waves, at intermediate times 
it may or may not approach a higher n > 0 stationary 
solution and remain for a quite macroscopic time in its 















FIG. 4: Snapshots of the condensate density (thick solid blue line) at different evolution times (indicated in the panels) for a 
configuration with v = 0.75, C 2 = 0.3 and X = 2 that eventually reaches the n = 0 stationary solution (thin solid black line). 
The vertical dashed lines mark the boundaries of the supersonic region. The complete time evolution is shown in Movie 1 
Upper row: evolution of the system when the initial linear instability grows in the direction of the n = 0. Central and lower 
row: evolution of the system when the initial linear instability grows in the opposite direction. The complete time evolution is 
shown in Movie 2 
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FIG. 5: Similar to Fig. |4]but now for a system with v = 0.75, C 2 = 0.5 and X = 5, where the unstable mode presents a non-zero 
real part of the frequency so it oscillates while growing. The complete evolution is shown in Movie 3 


vicinity. As such solutions are dynamically unstable, the verge to the stable n = 0 solution after another stage of 
system must eventually depart from it and finally con- soliton and wave emission. 
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FIG. 6: Analog of Fig. |4]but for a system with v = 0.75, C 2 = 0.6 and A = 10 so there are two unstable modes. Upper and 
central rows: evolution of the system when the linear unstable mode n = 1 grows towards the n = 1 solution (depicted in a 
thin black line in the upper row) and remains there for some time until it collapses and evolves towards the n — 0 ground state 
(thin black line in the central row). The complete evolution is displayed in Movie 4 Lower row: evolution of the system when 
the linear unstable mode n = 1 grows in the opposite direction, emitting some waves and solitons before reaching the n = 0 
solution. The complete evolution is shown in Movie 5 


When A„(u, C 2 ) < X < X^_f.i/ 2 ('^X 2 ), some useful in¬ 
formation about whether the system does or does not 
spend time in the vicinity of a higher n > 0 stationary 
solution can be obtained along the lines of the discussion 
around Eq. of the previous subsection: depending on 
the sign of the projection of the initial noise configura¬ 
tion onto the n dominant unstable mode, the system can 
smoothly evolve towards the n non-linear solution or it 
can start evolving in the opposite direction. In this latter 
case, more solitons and waves need emitting to conserve 
the total number of particles so the evolution is more 
chaotic. Both trends are presented in Fig. [6l where we 
consider a simple configuration with only two unstable 
(n = 0,1) modes. In the upper and central rows, we 
show the case where the system smoothly evolves to the 
n = 1 solution, remaining there for a long time until it 
collapses and then reaches the ground state n = 0. In the 
lower row, we show the evolution when the state grows in 
the opposite direction, where after emitting some waves 
and solitons, the system reaches the ground state in a 
much shorter time scale. In both cases, after a (possi¬ 
bly very long) transient, the simulation ends up in the 
lowest energy n = 0 stationary solution. The complete 


evolution of both situations is displayed in Movie 4 and 
in Movie 5, respectively. 

This dramatic dependence on the initial conditions is a 
simplest example of the potentially chaotic dynamics of 
the system for Xi < X when many unstable modes are 
present. As another remarkable example of non-trivial 
behaviour. Movie 6 shows a case when the system hap¬ 
pens to intercept some other nonlinear stationary solu¬ 
tion, so that it remains trapped in a sort of metastable 
state in its vicinity for a macroscopically long time and 
the n = 0 solution is not yet reached for very long times 
on the order oi t = 10^. Even though all n > 0 station¬ 
ary solutions are dynamically unstable and the system 
must eventually depart from them, we cannot exclude 
that there exist regions of the (A, C 2 , v) parameter space 
for which the time scale of the decay to the n = 0 sta¬ 
tionary state might take arbitrarily long. 

B. Continuous emission of solitons 

This subsection expands the discussion of the CES 
mechanism briefly mentioned in Secs. IIII A II and IIII 
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and reports some among the main results of this work, in 
the direction of understanding such an intriguing regime. 
In contrast to the Bogoliubov-Cherenkov emission of 
sound waves that takes place in the upstream direc¬ 
tion [l^ , the continuous emission of solitons occurs in the 
downstream direction towards x = -foo (Hill. Similar 
scenarios of soliton train emission have been predicted in 
Refs. [32| - [33 . and experimentally observed in In 
all these works the soliton emission process takes place 
in the vicinity of a localized defect potential where the 
condensate density is locally depleted and a localized su¬ 
personic flow appears. An interpretation of this physics 
in terms of a black-hole laser instability has been pro¬ 
posed in [ 3 ^. Related, but somehow different soliton 
emission mechanisms were discussed in [H, [s^. 

The first step in the direction of characterizing the CES 
regime is to identify in the (A, C 2 ,v) parameter space the 
regions of CES: cuts of such phase diagram along the 
(c 2 ,u) plane are shown in the different panels of Fig. [7] 
for growing values of X. The region where CES happens 
is indicated by black crosses and extends for growing X 
starting from the upper-left corner (i.e. the region with 
high V and low C 2 ). The CES regime is indeed qualita¬ 
tively related to the degree of instability of the central 
supersonic region, which is favored by a large flow speed 
V and a low C 2 . Even though a chaotic CES can be ob¬ 
served in very strongly unstable systems with many un¬ 
stable modes, we will not focus our attention on it and 
we will restrict the use of the CES expression to periodic 
soliton emission processes. 

As done in the previous section, to understand the 
physics underlying the CES process we start by consid¬ 
ering configurations with a single unstable mode where 
only the n = 0 nonlinear stationary solution is present. 
Several snapshots of the corresponding time-evolution are 
shown in Fig. [8] for a parameter choice that ends in the 
CES regime; the complete time evolution can be observed 
in Movie 7, In the upper row of the figure, we represent 
the initial evolution of the system: the emission of soli¬ 
tons right after the onset of the initial instability. In 
the lower row, we analyze in detail the CES mechanism, 
which persists in a periodic way indefinitely for arbitrar¬ 
ily long times. The soliton that tries to be emitted in 
the upstream direction is dragged by the flowing conden¬ 
sate and bounces back. Eventually, it ends up in the 
downstream region and travels towards x —>■ + 00 . After 
the soliton has gone, the density modulation in the su¬ 
personic region begins to grow again until a new soliton 
is generated. Continuous periodic repetition of this pro¬ 
cess leads to the emission of a train of solitons into the 
downstream region. 

In order to check quantitatively the periodicity of 
the soliton emission process, we introduce the quantity 
tsoi{N), which is defined as the time lapse between the 
emission of the two consecutive {N — l)th soliton and 
iVth soliton. If the emission of solitons is periodic, tsoi 
should be constant and equal to the period of emission T, 
tsoi = T. For counting solitons, we monitor the density 


at a generic point xq ^ A/2 in the downstream region as 
a function of time, as shown in the inset of Fig. [HI The 
minima of the density correspond to the passage of a soli¬ 
ton. In the main panel we represent the time-evolution 
of the time lapse tsoi as a function of the number of emit¬ 
ted solitons, Nsoi- After some irregular transient, we see 
that this quantity approaches a constant value T, mean¬ 
ing that the soliton emission process becomes an almost 
perfectly periodic one. 

Interesting light on the soliton emission process can 
be obtained by studying the dependence of the emission 
period T on the system parameters. As an example, in 
Fig. [ini we plot the CES period as a function of C 2 for 
different values of v for fixed X = 2. Every curve stops 
at some critical value of C 2 : after this critical value the 
system no longer performs a CES and a single soliton is 
in fact emitted to the upstream region while the system 
relaxes to the n = 0 nonlinear solution. Below this criti¬ 
cal value, the period increases with C 2 as expected since 
for larger values of C 2 the instability of the the supersonic 
region becomes weaker. This physical interpretation also 
explains the decreasing value of T for increasing v seen 
in Fig. 1101 this argument also agrees with the numerical 
observation of the decrease of T for increasing X. 

As a last ingredient before proposing a qualitative ex¬ 
planation of the CES mechanism, it is useful to remind 
a few crucial features of dark solitons in atomic conden¬ 
sates. As we have seen in the previous sections, dark soli¬ 
tons are often emitted into the upstream and/or down¬ 
stream external regions to compensate for the increased 
density in the internal region. In the external region, the 
coupling constant is homogeneous g{x) = 1 and the con¬ 
densate flows at a constant speed u, so dark solitons are 
described by solutions of the GP equation of the form 


«'(x) = e™" 


ik + \/l — tanh(-\/l — x) , 


(24) 


parametrized by a single real parameter — 1 < fc < 1. 
The density profile shows a dip with a minimum den¬ 
sity equal to k^ and, apart from some irrelevant global 
time-dependent phase, it rigidly moves at a fe-dependent 
velocity 


Vs = V -\- k (25) 

which results from the Galilean combination of an overall 
drag at the condensate speed v and a relative motion at 
speed k equal to the local value of the speed of sound at 
the density minimum: the deeper the soliton, the smaller 
its relative velocity with respect to the surrounding fluid. 

A soliton is typically emitted into the upstream direc¬ 
tion in order to conserve the total number of particles. 
We can expect that its minimum density k^ is a roughly 
decreasing function of the amplitude of the non-linear 
stationary solution n = 0. As the soliton moves in the 
upstream direction, the sign of k is fixed to fc < 0. In 
particular, we have numerically observed that the k pa¬ 
rameter of the emitted soliton depends very weakly on 
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FIG. 7: Phase diagram in the {c 2 ,v) plane for different values of X. The blue dots represents numerical simulations that end 
in the n = 0 ground state (labeled as GS in the legend). The black crosses represents simulations in which CES has been 
observed. The solid curves from bottom to top represent the lower boundary of the linear instability regions defined in Eq. [3 
with n = 0 (blue), 1/2 (green), 1 (red), 3/2 (pink), etc. The dashed oblique straight line is the upper boundary of the region 
where the flow is everywhere subsonic. 


the initial condition so we can consider it to be a func¬ 
tion of the system parameters, k = k(v,C 2 ,X). As the k 
parameter fixes the relative velocity of the soliton with 
respect to the underlying fluid, a qualitative change of 
behavior is expected to occur around v = \k(v,C 2 ,X)\. 
For V < |fc(?;,C 2 ,Al)|, the soliton can freely travel in 
the upstream direction. However, in the opposite case, 
V > |fe(u, C 2 ,^)|, the flow speed v is larger than the rel¬ 
ative soliton velocity |fc|, so the soliton cannot travel to 
the upstream region and is dragged by the condensate 
flow into the downstream direction, originating the CES 
mechanism previously described. 

More insight on the transition to the CES regime 
can be obtained through Movies 8, 9 and 10. These 
movies report the result of numerical simulations for 


fixed V = 0.8, C 2 = 0.4 and growing X = 2.2 (Movie 8), 
X = 2.4 (Movie 9), and X = 2.5 (Movie 10). As X is 
increased, the soliton must become deeper and deeper in 
order to allow the system to relax to the n = 0 nonlinear 
stationary state, which implies a reduction oi k. As a 
result, the soliton becomes slower and slower: while for 
X = 2.2 and A = 2.4 the dark soliton is able to form 
and escape to the upstream direction so to compensate 
for the higher density of the n = 0 nonlinear stationary 
solutions, this can no longer happen for X = 2.5, where 
we have in fact entered the v > k regime and the soli¬ 
ton is forced to bounce back to the central region. As 
a result, it finally escapes in the downstream direction 
leaving the system again in an unstable configuration, 
so that the CES process keeps periodically repeating for 
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FIG. 8: Snapshots of the condensate density (solid blue line) at different evolution times (indicated in the panels) for a 
configuration with v = 0.95, C 2 = 0.4 and X = 2 that eventually reaches the CES regime. The horizontal black line represent 
the initial homogeneous condensate density. After some transient, the continuous emission of solitons begins. Upper row: initial 
evolution of the system when the initial linear instability sets in. The vertical black line marks the x = 0 position. Lower row: 
periodic evolution at late times after the system has reached the CES regime. The vertical dashed lines marks the boundaries 
of the supersonic region. At t = 1940, a soliton emerges near the black hole (subsonic-supersonic) interface at a: < 0. This 
soliton cannot travel upstream and then it bounces back, traveling to the downstream region and finally crossing the supersonic 
region, see plot at t = 1950. After this soliton has left around t = 1968, the density grows again and a new soliton starts 
growing near the interface. Finally, at t = 1982, the system recovers the same state as at t = 1940 and the process repeats. 
The complete time evolution is shown in Movie 7 



Nsol 


FIG. 9: Main panel: plot of tsoi, the time lapse between 
the arrival of consecutive solitons, as a function of the soli¬ 
ton number: after the emission of some solitons, the system 
reaches a regime of periodic soliton emission of period. Inset: 
time evolution of the density at a given point, p{xo,t) with 
*0 = 50. The density minima correspond to the periodic pas¬ 
sage of solitons. System parameters: v = 0.95, C 2 = 0.4 and 
X = 2. 

indefinite times. 

While this qualitative observation is in good agreement 



FIG. 10: Soliton emission period T as a function of C 2 for 
(from top to bottom) v = 0.8, 0.85, 0.9, 0.95 and for fixed X — 

2 . 


with the numerical results, it is far from offering a com¬ 
plete theoretical picture. In particular, finding a quan¬ 
titative relation between the minimum amplitude of the 
soliton emitted upstream and the parameters of the sys¬ 
tem, i.e., the function k{v,C 2 ,X), is still a theoretically 
unsettled problem. Obtaining and understanding a re¬ 
lation of this kind is an open question that should be 
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addressed in future works. 

The CES behavior in the X > Xi(u, C 2 ) regime with 
several unstable modes can be much more complex due 
to the presence of other n > 0 nonlinear stationary solu¬ 
tions. Typically, the system emits several solitons before 
reaching the periodic CES regime described previously. 
We show this scenario in Movie 11, which summarizes 
the result of a simulation for v = 0.65, C 2 = 0.3,X = 8. 
In other situations, the system first evolves towards a 
n > 0 solution, remains in the vicinity of such solution 
for quite some time, then departs again from it and, after 
some transient, reaches the CES regime. We show an ex¬ 
ample of this latter case in Movie 12: in this simulation, 
the system stays in the vicinity of the dynamically unsta¬ 
ble n = 1 solution for quite some time before entering in 
the CES regime. Finally, for a large number of unstable 
modes, the periodicity of the soliton emission may com¬ 
pletely disappear, with the system strongly oscillating 
around different solutions with n > 0. This chaotic sce¬ 
nario is illustrated in Movie 13, where the system shows 
a very complicated aperiodic emission of solitons. 


IV. DISCUSSION AND COMPARISON WITH 
OPTICAL LASERS 

On the basis of the results presented in the previous 
Sections, we are now in a position to critically discuss 
the actual physical content of the widely used expression 
hlack-hole laser. 

Historically, when this expression was first introduced 
by Corley and Jacobson [I|, the authors had in mind the 
linear dynamical instability of the configuration with a 
pair of neighboring black- and white-hole horizons. In 
terms of laser devices in quantum optics, this dynami¬ 
cal instability corresponds to the linear instability of the 
electromagnetic vacuum state of the cavity when the (un¬ 
saturated) gain exceeds losses. In this case, any weak sig¬ 
nal due to, e.g., spontaneous emission, is able to trigger 
the instability and quickly gets coherently amplified by 
stimulated emission. 

The actual operation of a continuous-wave laser de¬ 
vice [H takes place in a very different regime, where 
nonlinear gain saturation brings the system to a non¬ 
linear stationary state where (saturated) gain and losses 
exactly compensate each other and the field keeps os¬ 
cillating at a well-defined frequency with a well-defined 
amplitude. This clean and coherent laser emission is at¬ 
tained at long times after switch-on once all (possibly 
complex) transient phenomena are gone: from the point 
of view of the field dynamics, the long-time limit is a limit 
cycle in the classical dynamics, and the only observable 
quantum effect are small fluctuations in the oscillation 
frequency, leading to a very long but finite coherence 
time of the emission: in the simplest cases, this slow 
decoherence of the laser emission goes under the name of 
Schawlow-Townes linewidth. 

Even though from a quantum optics perspective the 


expression black-hole laser suggests a spontaneously os¬ 
cillating behavior with no memory of the initial state 
and strongly determined by nonlinear effects, the com¬ 
plex nonlinear dynamics of analog models in the two- 
horizon configurations has began being investigated only 
very recently in Building on top of these works, 

we have seen in the previous sections that the long-time 
behavior of such configurations is much richer than the 
one of standard optical lasers. 

In many cases (Sec. IIII All , the system is in fact able 
to quickly find a new time-independent stationary state 
where the instability is quenched by a suitable redistri¬ 
bution of the atomic density which eliminates the super¬ 
sonic character of the central region. During this rapid 
“evaporation” of the horizons, the extra energy and par¬ 
ticle density are compensated by the emission of a short 
sequence of sound and soliton pulses in the condensate, 
whose properties strongly depend on the (classical or 
quantum) noise triggering the initial linear instability. 

On the other hand, regimes fSec. IIII Bl) where the den¬ 
sity and the current keep oscillating in a periodic way in 
time can be considered as a hydrodynamic counterpart 
of a continuous-wave monochromatic optical laser oscil¬ 
lation, leading to a continuous and periodic emission of 
sound waves in the condensate in the form of a train of 
solitons. The excellent periodicity and insensitivity to 
initial conditions of this soliton train guarantees a very 
high degree of coherence of the emission, closely analo¬ 
gous to the one of the output beam of an optical laser 
device. On this basis, a rigorous reader may then argue 
that the black-hole lasing expression should only be used 
to refer to this latter case. 

As a final remark, it is worth noting that in contrast 
to standard single-mode laser devices which emit light 
of a single frequency, the solitons emitted by the black- 
hole laser device are intrinsically nonlinear objects con¬ 
taining a number of frequency components: differently 
from pulsed lasers where the pulse-generating mechanism 
serves to modulate a much faster carrier frequency and 
creates tightly spaced sidebands, here the nonlinearity 
provides new frequencies at integer multiples of the fun¬ 
damental one. In a qualitative picture, these components 
can be understood as resulting from a phase-locking of 
the different modes of the central, super-sonic region that 
are simultaneously oscillating. The energy for these os¬ 
cillations is of course provided by the macroscopic flow of 
the condensate, and the mechanism responsible for the 
phase-locking of the different components originates from 
the hydrodynamic nonlinearities at large density modu¬ 
lations. 


V. CONCLUSIONS AND OUTLOOK 

In this work we have reported the results of an ex¬ 
tensive campaign of numerical simulations of the Gross- 
Pitaevskii equation describing the time-evolution of an 
atomic Bose-Einstein condensate. For the sake of con- 
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ceptual simplicity, our numerical study was performed 
on a theoretically most convenient two-horizon configu¬ 
ration, which is theoretically convenient since it allows 
to suppress additional phenomena such as Bogoliubov- 
Cherenkov emission of sound waves. While this config¬ 
uration may require a challenging experimental effort to 
simultaneously control in a spatially selective way both 
the atomic interaction constant and the external poten¬ 
tial, we expect that the qualitative conclusions of our 
work extend to a wider range of black-hole laser configu¬ 
rations, e.g. the ones used for the experimental observa¬ 
tion of the (linear) black hole instability in [T^. Prelim¬ 
inary numerical work for this specific configuration 
suggests that, at least at short times, the black-hole las¬ 
ing mechanism is taking place along very similar lines 
also in this case. 

In agreement with existing theoretical and numerical 
results about linear instabilities and nonlinear stationary 
solutions [ 21 |,[ 23 , our simulations provide useful physical 
insight on the system behavior as a function of the initial 
system parameters. In some cases, the central unstable 
region is eventually evaporated away and the system at 
late times approaches a stationary stable configuration. 
In other cases, the system approaches a novel regime of 
continuous and periodic emission of solitons for indefinite 
times. Interesting analogies and differences between the 
continuous soliton emission in black-hole lasers and the 
operation of an optical laser device can be drawn from 
our simulations. On this basis, a more restrictive and co¬ 
herent use of the expression hlack-hole laser is proposed. 

As a further remarkable result, our simulations high¬ 
light how the system typically forgets its initial state and, 
after possibly very different transients, the late time dy¬ 
namics is only determined by the parameters of the con¬ 
figuration and not by the initial condition. This fact pro¬ 
vides numerical evidence of the existence of a no-hair the¬ 
orem also for analog black-hole lasers, in a similar fashion 
to the no-hair theorem for analog black holes [H, . 

While all our discussion has focussed on the classical 
black-hole laser dynamics, future work will address the 
effect of quantum fluctuations on these systems: taking 
again inspiration from quantum optics and laser theory, 
we may legitimately expect that quantum effects should 
reduce to weak fluctuations in the emission frequency 
and, therefore, to a slow decoherence of the soliton emis¬ 
sion analogous to the long, but finite coherence time of 
optical laser devices [Hi. 

From the experimental point of view, we can reason¬ 
ably expect that a n up graded version of the ultracold 
atom experiment in |15| with a longer condensate sample 
will soon be able to investigate the nonlinear dynamics 
of the black-hole laser at late times after the onset of the 
dynamical instability and hopefully characterize the rich 
phenomenology discussed in this work. 

Even though all the discussion in this work has been 
carried out having an atomic implementation in mind, of 


course we expect that most conclusions can be directly 
transferred to analog models based on fluids of light in 
the so-called propagating geometry, for which theoreti¬ 
cal works on black hole configurations have recently ap¬ 
peared [ 4 II, 113 as well as first experimental evidences of 
superfluid behaviors (d^ . 

Apart from analog gravitational purposes, the black- 
hole laser setup here described and, in particular, the 
characterized non-linear soliton laser regime, can provide 
an interesting scenario for the investigation of quantum 
transport, within the emergent field of atomtronics 
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Appendix A: Numerical scheme 


The GP equation, Eq. (EJ, is integrated making use of 
a time-splitting finite-difference (TSFD) scheme pa . \4l\ . 
For this purpose, we write the GP equation as 


' dt 


HGp{x,t)'i/{x,t) 

[Ho{x,t) + Hc{x,t)]'i>{x,t) (Al) 


where Hq is the kinetic term and Hq = g{x){\'^{x,t)\‘^ — 
1) is the term associated with the contact interaction. We 
divide the spatial interval [—Lg/ 2 , Lg/ 2 ] in steps of size 
Ax and impose periodic boundary conditions. The time 
interval [0, t] is divided in n steps of size At . Within this 
scheme, the time evolution operator for time tk = kAt 
to tfc+i = (k + l)At can be approximated as: 




(A2) 


where^ = g{x){\^{x,tk)\'^ - 1), ^ix,tk) = 

e~'^^'^^'i'{x,tk)- The total evolution operator can be 
then written in the form: 
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= U{t)'i'{x,0) 
U{t) = 


In the TSFD method, instead of using the usual Fourier 
transform to compute the kinetic evolution operator 
^-iHoAt^ we use the Crank-Nicolson method within a fi¬ 
nite difference scheme The main advantage of 

using the Crank-Nicolson method is that it is uncondi¬ 
tionally unstable, which allows us to use higher values of 
Ax, At. 

In addition, such a method keeps working even when 
the kinetic term involves derivatives with a spatially de¬ 
pendent coefficient. This feature will be crucial to imple¬ 
ment the absorbing mechanism needed to suppress arti¬ 
facts in long-time simulations: our strategy to this pur¬ 
pose is to add a diffusive term on the boundaries of the 
integration box so to absorb all emitted perturbations 
and avoid their return to the central region of interest. 
Specifically, we replace Hqp by Hp = Hqp + iHa in Eq. 
(EB), where the absorbing term Ha has the form: 

HA{x,t) = G{x) - Eo[|vI/(x, 1)1^ - 1]) (A4) 

in terms of a function G{x) localized on the boundary of 
the grid and sufficiently smooth in order to avoid reflec¬ 
tions. In particular, this term vanishes when acting on 
the unperturbed homogeneous plane wave . Within 
Ha, the first term between the brackets of the r.h.s is 
a diffusive term that absorbs the Fourier components 
outside the plane wave while the second term is a 
source term that serves to compensate the undesired loss 
of atoms introduced by the diffusive term. 

The new diffusive GP equation can still be integrated 
using the mentioned TSFD method by just replacing Hq 
and Hq in Eqs. (THIl-dMl) hy Hp) and Hx, where Hp 
is the term that contains the derivatives (the kinetic and 
diffusive terms) and Hx is the term that takes into ac¬ 
count the non-linearities (the interaction and the source 
term), i.e., 

Hx = [g{x) - iFoG{x)][\'i>{x,t)\‘^ - 1] 

As the source term does not conserves the norm, it 
is useful to refine the simple scheme discussed above 
by introducing a predictor-corrector method in order 
to provide a better estimation for the non-linear term 
in Hx, similar to that of Refs, [s^, [Ml. At every 
step, we perform a first iteration to compute 4''(x,tfc) = 
tfc) where '^{x,tk) is now '^{x,tk) = 

e~^^°^'S{x,tk) andH^ = [g{x)-iFoG{x)]{\^{x,tk)\‘^- 
1). After this first iteration, we replace ^{x,tk) by 


\k=0 


r 


(A3) 


'^'{x,tk) -f ^ix,tk) 


/2 in H^ and we perform a second 


iteration to obtain the final value of ^'(a;,tfc). 

Typically we take for G{x) a Gaussian shape with a 
typical width on the order of 10^ and values in the range 
I — 10 for Fq and Dq. With parameter choices in such in¬ 
tervals, we have not seen any significant variation in the 
result of the simulations. As a further check, we have 
verified that the diffusive scheme does not introduce any 
spurious result by comparing with the result of integrat¬ 
ing the plain GP equation without diffusion but on a 
much larger spatial grid. 

It is worth noting that such a refined treatment of the 
boundaries is essential for the study of the CES regime 
of Sec. IIII Bl Indeed, the observation of the perfect peri¬ 
odicity of the soliton emission is in fact sensitive to small 
spurious reflections at the boundaries. As in all other 
cases, such artifacts are of course further suppressed with 
a sufficiently wide integration box. 


Appendix B: Movies 

In this Appendix, we summarize the system parame¬ 
ters used for the Movies that are discussed in this work. 
In all the videos, the upper panel displays the time evolu¬ 
tion of the density profile (solid blue line) while the lower 
profile shows that of the sound (solid blue) and flow (red 
dashed) velocities. The limits of the internal region are 
marked by two vertical dashed black lines. 

• Movie 1: V = 0.75, C 2 = 0.3 and X = 2 satisfying 
A'o('c,C 2 ) < X < Xi/ 2 {v,C 2 ). The unstable mode 
initially grows in the same direction of the n = 0 
nonlinear stationary solution (solid black line in the 
upper panel). 

• Movie 2: Same parameters as Movie I but now the 
unstable mode initially grows in the opposite direc¬ 
tion. 

• Movie 3: v = 0.75, C 2 = 0.5, A = 5, with 
Ai/ 2 (u,C 2 ) < a < Ai(u, C 2 ). We observe the oscil¬ 
lating behavior during the growth of the unstable 
mode, eventually reaching the n = 0 ground state. 

• Movie 4: v = 0.75, C 2 = 0.6, A = 10, with 
Ai(u, C 2 ) < A < A 3 / 2 (u,C 2 ). The unstable mode 
grows towards the n = I stationary solution. Since 
that solution is dynamically unstable, the system 
stays in its vicinity for a large but finite time, then 
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it departs from it to finally reach the n = 0 ground 
state (solid black line in the upper panel). 

• Movie 5: Same parameters as Movie 4 but the ini¬ 
tial unstable mode grows now in the opposite di¬ 
rection. The system emits some waves and soli- 
tons and then reaches the n = 0 ground state in a 
shorter time. 

• Movie 6: v = 0.9, C 2 = 0.75, X = 20, with 
X 5 / 2 (u, C 2 ) < X < X 3 (u, C 2 ). After some transient, 
the system reaches a configuration near the n = 1 
stationary solution (black line in the upper panel), 
around which the system is still oscillating at the 
end of the simulation {t = 10000). 

• Movie?: v = 0.95, C 2 = 0.4, X = 2, with 
Xo(t’,C 2 ) < X < Xi/ 2 (n,C 2 ). Instead of evolving 
to the ground state n = 0 (solid black line in the 
upper panel), the system reaches a regime of con¬ 
tinuous emission of solitons. A slow-motion effect 
has been introduced at large times in order to make 
easier the observation of the mechanism of periodic 
emission of solitons. 

• Movie 8: v = 0.8, C 2 = 0.4, X = 2.2, with 
Xo{v,C 2 ) < X < Xi/ 2 (u, C 2 ). We observe that the 
system emits a soliton to the upstream region and 
reaches the ground state n = 0 solntion (solid black 
line in the upper panel). 

• Movie 9: Same parameters as Movie 8, but with 
X = 2.4. The soliton is still able to travel in the 


upstream direction. 

• Movie 10: Same parameters as Movie 8, but with 
X = 2.5. In this case, the upstream traveling soli¬ 
ton bounces back and the system reaches the CES 
regime. 

• Movie 11: V = 0.65, C 2 = 0.3, X = 8, with 

Xi(n, C 2 ) < X < X 3 / 2 (u, C 2 ). Even when more 
than one unstable mode is present, after some tran¬ 
sient, the system reaches the same CES regime as 
before. 

• Movie 12: V = 0.65, C 2 = 0.1, X = 8, with 

X 5 (u, C 2 ) < X < Xii/ 2 (u, C 2 ). The system evolves 

towards the n = 1 non-linear solution (black line in 
the upper panel). As this solution is dynamically 
unstable, after some time the system departs from 
it and eventually reaches the CES regime. 

• Movie 13: V = 0.9,C 2 = 0.2,X = 20, with 

X 5 (u, C 2 ) < X < Xii/ 2 {v,C 2 ). Due to the high 
number of unstable modes, the system is contin¬ 
uously emitting solitons but in a chaotic way, in 
contrast to the usual periodic CES regime. The 
non-linear n = 5 stationary solution is depicted as 
a solid black line in the upper panel. 
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